#include <stdio.h>
	#include <math.h>
	#include <string.h>
	#include "helper.h"
	#include "init.h"
	#include "sor.h"
	#include "uvp.h"
	#include "boundary_val.h"
	
	int main(int argn, char** args){
	
		/*
		 * - read the program configuration file using read_parameters()
		 */
	
		double Re, UI, VI, PI, GX, GY, t_end, xlength, ylength, dt, dx, dy, alpha,omg, tau,eps,
		dt_value,dp,y1,y2,dt_path,dt_streak,dt_insert;
		int ** Flag;                   
		int  imax, jmax, itermax,wl,wr,wt,wb,N,x1,x2;
		char szFileName[]="step.dat";
		int fluidCellCount=0;
		char * problem= (char*)malloc(20);
		read_parameters(szFileName,&Re,&UI, &VI, &PI,&GX,&GY,&t_end,&xlength,
				&ylength,&dt, &dx, &dy,&imax, &jmax,&alpha, &omg, &tau,&itermax, &eps,
				&dt_value,&wl,&wr,&wt,&wb,problem,&dp,&N,&x1,&y1,&x2,&y2,&dt_path,&dt_streak,&dt_insert);
	
		/*
		 *  - set up the matrices (arrays) needed using the matrix() command
		 */
	
		double **U, **V, **P,**F, **G, **RS;
		//Added
		double **Psi;
		double **Zeta;
		
		 struct particle_line *streakline;
	     struct particle_line *pathline;
		 int insrt;
		 const int partmax =1000;
		 
		 
		U =matrix(0, imax+1, 0, jmax+1);
		V =matrix(0, imax+1, 0, jmax+1);
		P =matrix(0, imax+1, 0, jmax+1);
		F =matrix(0, imax, 0, jmax);
		G =matrix(0, imax, 0, jmax);
		RS=matrix(0, imax+1, 0, jmax+1);
		/*Initialize the flag matrix*/
		Flag=imatrix(0,imax+1,0,jmax+1);
		//ADDED
		Psi = matrix(0,imax+1,0,jmax+1);
	   Zeta = matrix(0,imax+1,0,jmax+1);
		
		
		/*
		 *  - create the initial setup init_uvp(), init_flag(), output_uvp()
		 * init_flag() ??
		 */
	
		init_matrix(F,0,imax,0,jmax,0);
		init_matrix(G,0,imax,0,jmax,0);
		
		
		init_uvp(UI, VI, PI, imax, jmax, U, V, P,problem);
		
		init_flag(problem,imax,jmax,Flag);
	//	pathline=set_particle(N, x1, y1, x2, y2);
		//streakline=set_particle(N, x1, y1, x2, y2);
		char *szProblem=(char*)malloc((strlen(szFileName)+1)*sizeof(char));
		strcpy(szProblem, szFileName);
	    //Calculate the fluid cells
		fluidCellCount=getFluidCellCount(imax,jmax,Flag);
	
		/*
		 * - perform the main loop
		 */
		double t=0, res;
		int n=0, it=0, counter=1;
		while (t<t_end) {
	
			//calculate_dt(Re, tau, &dt, dx, dy, imax, jmax, U, V);
			//pathlines(pathline, N,dt_path, t, dx, dy,dt, imax, jmax, U, V,Flag,szProblem,counter);
			
			
			//pathlines(pathline,N,szProblem,dt_path,t,dx,dy,dt,imax,jmax,U,V,Flag);
           // streaklines(streakline, N,dt_insert, dt_streak, t,dx, dy, dt,imax, jmax, U,V,Flag,szProblem,counter);
			dt = 0.003;
			boundaryvalues(imax, jmax, U, V,wl,wr,wt,wb,Flag);		
			calculate_fg(Re, GX, GY, alpha, dt, dx, dy, imax, jmax, U, V, F, G,Flag);
			calculate_rs(dt, dx, dy, imax, jmax, F, G, RS,Flag);
			res=1.0;
			it=0;
			while (it<itermax && res>eps) {
	
				sor(omg, dx, dy, imax, jmax, P, RS, &res,Flag,fluidCellCount,dp);
				
				it++;
			}
	
			//printf("SOR iteration count: %d\n", it);
	
			calculate_uv(dt, dx, dy, imax, jmax, U, V, F, G, P, Flag);
			//ADDED
		//	calculate_psi_zeta(dx,dy,imax,jmax,U,V,Psi,Zeta,Flag);		
			// Write it to the file in every 100. iteration and t += 0.5 seconds
			if ( (n%18) == 0){
		//	if (t>=t_end-0.1*dt_value){
				printf("Writing to file: Iteration # %d, output filename ends with # %d, at time t=%2.3f\n"
				,n, counter, t);
				write_vtkFile(szProblem, counter, xlength, ylength, imax, jmax, dx, dy, U, V, P,Flag);
				counter++;
				if (counter >100 ) return -1;
			}
	
			if (n ==100){
				printf("U[imax/2][7*jmax/8]...%f", U[imax/2][7*jmax/8]);
				
			}
			t=t+dt;
			n++;
		}
		
		printf("U[imax/2][7*jmax/8]...%f\n", U[imax/2][7*jmax/8]);
		
	
		free_matrix(U, 0, imax+1, 0, jmax+1);
		free_matrix(V, 0, imax+1, 0, jmax+1);
		free_matrix(P, 0, imax+1, 0, jmax+1);
		free_matrix(F, 0, imax+1, 0, jmax+1);
		free_matrix(G, 0, imax+1, 0, jmax+1);
		free_matrix(RS, 0, imax+1, 0, jmax+1);
	
	
		return -1;
	}
	